Predicting the dispersal and invasion dynamics of ambrosia beetles through demographic reconstruction and process-explicit modeling

Evaluating potential routes of invasion of pathogens and vectors of sanitary importance is essential for planning and decision-making at multiple scales. An effective tool are process-explicit models that allow coupling environmental, demographic and dispersal information to evaluate population growth and range dynamics as a function of the abiotic conditions in a region. In this work we simulate multiple dispersal/invasion routes in Mexico that could be taken by ambrosia beetles and a specific symbiont, Harringtonia lauricola, responsible for a severe epiphytic of Lauraceae in North America. We used Xyleborus bispinatus Eichhoff 1868 as a study subject and estimated its demography in the laboratory in a temperature gradient (17, 20, 26, 29, 35 °C), which we then used to parameterize a process-based model to estimate its metapopulation dynamics. The maximum intrinsic growth rate of X. bispinatus is 0.13 with a thermal optimum of 26.2 °C. The models suggest important regions for the establishment and dispersal the states of Veracruz, Chiapas and Oaxaca (high host and secondary vectors diversity), the Isthmus of Tehuantepec (connectivity region), and Michoacán and Jalisco (important avocado plantations). The use of hybrid process-based models is a promising tool to refine the predictions applied to the study of biological invasions and species distributions.

www.nature.com/scientificreports/COMPADRE [https:// compa dre-db.org/ Data] only includes two species of Curculionidae, S. ventralis and H. hampei) or are not entirely freely accessible.Regarding the accessibility of information where software does exist 105 , it may only support a subset of methods, and a lack of technical knowledge can impede their use, thus limiting potential applications.However, there is optimism due to the increase in data collection in databases in recent years.Coupled with the availability of freely accessible environmental layers at fine resolutions, current computing power, and growing knowledge applied to software development, these factors are expected to contribute significantly to overcoming the aforementioned shortcomings.
Although explicit examples that leverage morphological, phylogenetic, or environmental information among species to categorize them as analogous are few, there is implicit evidence in many works indicating the closeness they share in several of these aspects.An instance of this phenomenon is observed with ambrosia beetles.In recent years, it has been confirmed that the sister species Xyleborus affinis Eichhoff 1868, X. volvulus (F.) 1775, and X. bispinatus Eichhoff 1868 (until 2006 not differentiated from X. ferrugineus 106 ) can function as effective secondary vectors of H. lauricola 79,98,[107][108][109][110][111][112] .Among the three species of ambrosial beetles, X. bispinatus stands out due to its extensive geographic distribution across the American continent and recent occurrence in Europe and New Guinea 113 , its ability to survive in a diverse range of hosts (including avocado) 107 , and its association with the aforementioned phytopathogenic fungus 114 .Furthermore, considering its phylogenetic proximity 102 and environmental affinity 25,115 with X. glabratus, we propose that this species could serve as a model for estimating the demography of X. glabratus and suggesting management plans, particularly given the phytosanitary constraints associated with working with this exotic species.
Eradication of an invasive species over large areas is extremely difficult when it has been established for a long time 43,116,117 .The alternative approach is to predict and understand invasion pathways and processes in order to prioritize strategies to control its arrival and spread 118,119 .In this study, we use an explicit process-based dynamic model to estimate the processes and resulting geographic patterns of ambrosia beetles in Mexico.We parameterize the model using demographic rates obtained from laboratory rearing for X. bispinatus and niche modeling.We demonstrate how simulations from process-based models can complement phytosanitary tools by considering environmental, demographic, and dispersal aspects.This approach proves valuable for suggesting monitoring and early detection protocols to contribute to decision-making in programs for the prevention, management, and eradication of invasive species.

Laboratory
Adjustment of the growth function and optimal rate of growth Survivorship curves and oviposition rates of X. bispinatus were obtained to estimate the intrinsic growth rate r as function of temperature and the optimum temperature t opt at which the maximum growth rate (r max ) occurs, and kurtosis of the fitted curve (v 0 ) which we interpreted as a measure of niche breadth.r max occurs at 26.26 °C t opt .r max = 0.13 individuals/individuals/day (Table 1).The coefficient of determination of the curve was R 2 = 0.67, while the hypothesis tests for the non-linear regression parameters showed a high level of significance (Fig. 1; Table 1).www.nature.com/scientificreports/

Dispersal dynamics modeling
Environmental suitability and correlative model using the intrinsic growth rate r The selected correlative model incorporated environmental layers bio4, bio9, bio10, bio12, bio15, bio17, and bio18 120 and satisfied the evaluation criteria based on the following parameters: omission rate in test data of 0.074, omission rate in training data of 0.063, average omission rate of training and test data of 0.068, background prevalence of 0.816, and passed the binomial test (P-value = 0.029).This model showed sites with environmental suitability for the species along the slope of the Gulf of Mexico and in some states on the Pacific coast (Fig. 2A).The r i map (Fig. 2B) showed areas with different levels of suitability in the states of Veracruz, Oaxaca, Tabasco, and Chiapas, especially important due to the occurrence of higher numbers of species of Lauraceae hosts and secondary vectors, as well as in the Isthmus of Tehuantepec, which also connects the western and eastern regions of the country.Highly suitable climates were also recorded in the Yucatan peninsula and the avocado-producing area comprised of Michoacán, Jalisco, Nayarit, part of the State of Mexico, Morelos, and Colima.Towards the north of Mexico, a decrease in r i values were observed, which could mean a natural barrier for the dispersal and establishment of the species at these latitudes (Fig. 2B).

Metapopulation process-explicit model
The states of Veracruz, Chiapas, Oaxaca, and the Isthmus of Tehuantepec (including Tabasco) and the central avocado region of the country were the most threatened areas requiring monitoring for known or secondary vectors and the disease (also see synthesis of the simulations dynamics in Supplementary Table S2).Different spread speeds were observed in the invasion dynamics depending on the starting points of the simulation, the fastest being those that began from the ports of Salina Cruz and Veracruz, which took approximately 650 days to cover these areas (Fig. 3), followed by the ports of Manzanillo and Altamira with 760 and 970 days, respectively (Fig. 4).In Manzanillo, the invasion of the avocado area was the most important and rapid, completed within 350 days.In the case of Nogales customs, the elapsed time to cover these areas was 2670 days, while from Texas (a site where X. glabratus and X. bispinatus coexist and H. lauricola is known to be present) was 1700 days (Fig. 5).
Finally, the simulation from the port of Ensenada in Baja California did not prosper due to a lack of suitable sites (see all simulations in Supplementary Videos online).

Optimal growth temperatures, maximum fitness and growth function
Insects are susceptible to thermal variation [121][122][123][124][125][126][127][128][129] .In artificially reared Scolytinae beetles, temperature has been documented to be a key factor for their physiology, abundance, and distribution [130][131][132][133][134][135][136][137][138] as well as that of their fungal symbionts 67,139,140 .Here, rearing of X. bispinatus yielded similar results in terms of demographic www.nature.com/scientificreports/parameters relative to its artificially reared relatives at comparable temperatures 132,134,[141][142][143][144][145] .Our experiment showed the highest population growth between 26 and 29 °C.This optimal range was also recorded in X. affinis, another native secondary vector that was reared following a similar methodology 136 and had substantial population growth between 20 and 30 °C146,147 .In X. glabratus, Brar et al. 132 reported that the species completed its life cycle successfully at 24, 28 and 32 °C, with the highest rate of oviposition and development at 28 °C, and suggest through linear models that the lowest threshold for eggs and pupae would be at 13.8 °C and 11.1 °C, respectively.Likewise, these authors 132 reported that the rate of development observed for X. glabratus is similar to that reported in the literature for X. fornicatus and other Scolytinae species.
Although best growth for X. bispinatus occurred at 26 and 29 °C in the laboratory, the adjustment of the growth curve indicates an optimal temperature t opt at 26.26 °C, which might imply optimal breeding temperatures between 26 and 27 °C.Such values, however, were derived from individuals obtained from a population of Los Tuxtlas, Veracruz (Fadda et al., in review), and we should assume conservatism in this abiotic dimension of their fundamental ecological niche for extrapolation 25,115,[148][149][150][151] ; this is reinforced by data on the thermal niche of phylogenetically related species [152][153][154][155][156] .This implies that it would be possible to use demographic data of phylogenetically and functionally close species of Xyleborus as surrogates of their growth rates when there are phytosanitary restrictions for their direct use in the field or in laboratories that do not have the minimum necessary security levels 104,157 .

Mapping of intrinsic growth rate and influence of internal structure of the niche
Our methodological proposal suggests the use of demographic information from which optimal growth rates can be derived as a function of an abiotic factor (temperature in this case) and its expression in the geography (r i map).The r i map represents the relationship between the growth rate and the internal structure of the niche projected in geographic space 11,17 , thus indicating the distribution of hierarchies of climatic importance for the species where the process-explicit model will run.The r i map showed that the greatest demographic growth coincides with sites of greatest diversity of Lauraceae species 40 , including the municipalities of avocado production 94 (SIAP; http:// infos iap.siap.gob.mx/ gobmx/ datos Abier tos.php).This is a consistent result with what Lira-Noriega et al. 63 found using correlative ecological niche models for Lauraceae and Xyleborus species, with the neotropical region being the most vulnerable to the invasion of X. glabratus and Euwallacea spp., and for the distribution of X. bispinatus.Therefore, considering and incorporating the suitability expressed in the r i map as a function of one (e.g., temperature) or more dimensions of the ecological niche provides useful and pragmatic information that could be considered in decision-making for managing harmful species.Characterization the intrinsic growth rate is fundamental for forecasting future population trends 32 .www.nature.com/scientificreports/There are few antecedents about the use of species' growth thresholds and physiology-based correlative models in Xyleborus 158 .Most of this approach has focused on beetles of the genera Dendroctonus, Hylobius, Ips, and Hypothenemus 68 , which utilize phenological models integrating degree days to estimate the number of annual generations and predict population outbreaks of beetles.However, these models do not incorporate the upper and lower growth thresholds in these species, which poses a limitation when making predictions under changing climate scenarios 159 .An alternative approach that predicts this effect and uses correlative models is exemplified by the proposal of DeRose et al. 71 , who estimated the impact of increasing current and future temperatures on the distribution of D. rufipennis with BIOMOD.Evangelista et al. 61 used Maxent and Duehl et al. 69 CART (classification and regression tree) in D. ponderosae and D. frontalis, respectively, to estimate the suitability and importance of factors determining their distribution in North America.Jaramillo et al. 70 incorporated physiological data and degree days using CLIMEX to estimate the consequences of climate change on the number of generations of H. hampei (coffee berry borer) worldwide.These studies suggest that implementating processexplicit modeling represents a significant advancement in studying patterns and processes that influence the distribution of organisms with similar characteristics and of phytosanitary relevance.

Metapopulation process-explicit model
The implemented metapopulation process-explicit model, sometimes also referred to as hybrid models or coupled niche-population models 48 , combines information derived from correlative models with demographic and dispersal parameters to predict abundances and trajectories through a deterministic process.This model enabled us to test potential dispersal routes from multiple starting points across a large geographic extent.Although the predictions of our process-explicit simulations largely coincide with the suitabilities predicted by the correlative model, despite being initiated at different ports and customs across Mexico, they are useful for explaining likely patterns of movement across the territory of the ambrosia beetles (see Supplementary Videos online).The highest risk occurs if the simulation starts from Salina Cruz and Veracruz, probably due to the high environmental suitability and geographic proximity to the Isthmus of Tehuantepec and the avocado area and decreases in regions of northern Mexico with less suitability and connectivity (Fig. 3).Despite simulations initiated at Altamira and Manzanillo showing a slowdown in dispersion, their proximity to important avocado planting areas warrants monitoring (Fig. 4; Supplementary Table S2 online).The simulations from Nogales and from Texas take longer to arrive at highly suitable areas given the low intrinsic growth rates.However, from these sites the ambrosia beetles would be able to reach suitable sites such as Tamaulipas in the Gulf of Mexico or Jalisco in the Pacific (Fig. 5) from where a large expansion could occur in the interior of Mexico.Finally, the simulation from Ensenada did not prosper due to the low suitability, thus limiting the invasion progress to other regions.
The use of actual or hypothetical dispersal rates and of sites chosen for initiating the process-explicit simulations (e.g., ports, customs, or natural or human-facilitated dispersal rates) are helpful in assessing different invasion routes and the speed at which they can occur.This highlights the importance of connectivity corridors (Supplementary Video online).Additionally, we were able to identify source and sink sites/populations, which are often difficult to identify without estimations of population growth rates and environmental suitability throughout their distribution 13,30,160,161 .For our specific case, we assumed a constant dispersal rate of 5 km per day and a carrying capacity of 130,000 individuals per cell as hypothetical density dependency.However, such parameters could be adjusted based on available data from field or laboratory experiments or from rates reported in the literature.For example, Koch et al. 162 used a dispersal rate for X. glabratus in the USA of 54.8 km/year (approximately 0.15 km/day), which they estimated as the linear distance between infected counties between two years.Even though the dispersal rate estimated in this paper is hypothetical, determining the true dispersal of these ambrosia beetles in the field is difficult, and there is currently no information available regarding this aspect for X. bispinatus 163 .However, process-explicit models should allow for the evaluation of multiple dispersal rates/scenarios, perhaps due to the influence of anthropogenic interference and facilitation 53,63 .

Limitations and new horizons for the modeling of species niches and distributions
This study represents an effort to develop a tool based on ecological niche theory and process-explicit models that contribute to understanding the mechanisms and patterns responsible for the distribution of organisms.Despite illustrating its use in invasive species, its application can be extended to the study of other kinds of organisms and species assemblages at multiple scales in space and time 53 .The key advantage of process-explicit models lies in the capability of explicitly incorporating ecological mechanisms to better understand the determinants of species distributions 53,164 .However, the development of a process-explicit model depends primarily on the biological assumptions that are being explored and the mathematical form of the model.Some limitations of metapopulation models, such as the one implemented in this work, may include the mathematical model not being the most appropriate, challenges in obtaining demographic parameters, dispersion rates, and densitydependency factors, or inherent errors in the correlative models 54 used as input for estimating suitability and the map of the intrinsic rate of growth (r i ) 48 .
Most of the studies on biological invasions show the use of reaction-diffusion models that assume that space is continuous 165 .In contrast, our model uses discrete space, allowing for consideration of heterogeneous environments and the presence of dispersal barriers.In our approach, this can be addressed by assigning different dispersal rates (e.g., setting the dispersal rate between patch i and patch h to zero; δ ih = 0 ).A relevant factor not included in our study is the Allee effect, which plays a significant role as it emphasizes how regions of low suitability can act as biological barriers to dispersion 19,165,166 , and influence the temporal trajectory of invasion spread or population expansion.Moreover, Osorio-Olvera et al. 19 demonstrated that strong Allee effects can obscure the relationship between population abundance and niche structure because migration fails to start a population below an Allee threshold.However, this effect is less pronounced in our case because the www.nature.com/scientificreports/reproductive system of X. bispinatus is not affected by causes of inbreeding, and there is a low energy cost for finding reproductive partners [167][168][169] .Therefore, while the Allee effect may not be as relevant in the specific case of Xyleborus species due to their capacity for sib-mating or arrhenotokous inbreeding 67 , it is crucial to recognize its potential impact on other taxonomic groups that rely on it whenever this methodology and simulations are employed.Given the challenges associated with parameterizing demography-dependent process-explicit models, an alternative could be cellular automata models, where system states are presence/absence and a binarized niche model along with a first-neighbor dispersal kernel would suffice to estimate invasion paths 164 .However, these models cannot capture the level of complexity that process-explicit models like ours do, which incorporate demography and other key factors in many ecological and evolutionary processes.
We acknowledge that obtaining demographic data can be complex, especially for cryptic species that are challenging to study due to their complex behavior or lifestyle.Additionally, it may not always be feasible to follow cohorts.One way to overcome such limitation is through laboratory measurements where environmental conditions are controlled, although these come to the expense of not always reflecting the reality of the species' environment and the ways in which it affects its demography.Another limitation is related to the influence of several environmental dimensions on large-scale interpopulation variations.This might imply that demographic parameters estimated from one population may not be representative of the entire species 33 .Furthermore, the weighting of a few locally adapted units may underestimate the species' ecological niche 3 , especially if we lack information on whether these are source or sink populations 13,30,51,160 .However, laboratory studies and experiments are currently the most widely used methods in studying Xyleborus 108,130,134,136 and are considered a superior method for the study of ambrosia beetles 132 .This gains more importance considering that laboratory studies may be the only alternative for obtaining demographic information compared to other tools, such as COMPADRE and COMADRE (https:// compa dre-db.org/ Data), which may lack data for certain taxa and areas where they are distributed.
With regard to the use of correlative models of ecological niche, errors or biases can arise from various sources, including the quality of presence records (e.g., spatial, taxonomic, etc.), overall model quality, and the predicted suitability.Similarly, the scale and extrapolations in time and space can lead to over-or underestimations in these predictions 48,170,171 .We consider that these potential weaknesses could be mitigated because the training and validation of our model was carried out on a species within its native range under current climatic conditions.We also took careful consideration of errors on the input data and thoroughly evaluated the model.Furthermore, there is potential for improvement by incorporating microclimatic variables, which may be extremely important for the physiology of these species and their symbionts 67,[172][173][174] .Despite the potential drawbacks associated with the use of our process-based model, it remains a promising tool in ecology and evolution that could aid in better understanding species distributions.

Collection of individuals of X. bispinatus
Individuals were collected using various methodologies: infested trunks, females in flight, and bottle traps at the Estación de Biología Tropical of the Universidad Nacional Autónoma de México in the state of Veracruz.They were then transported to the Laboratorio de Entomología Molecular of the Instituto de Ecología, A.C., where they were conditioned and maintained in artificial culture media 141 at 26 °C and 60% relative humidity.The colony was maintained under these conditions until reaching the third filial generation in the laboratory, where live females were meticulously chosen for the experiment.

Culture medium
We used a modified culture medium of Biedermann et al. 141 : 45 g beech sawdust (Platanus mexicanus), 12 g agar, 6 g sucrose, 3 g casein, 3 g starch, 3 g yeast, 0.6 g Wesson's salt mixture, streptomycin 0.21 g, 1.5 mL of wheat germ oil, 3 mL of 96% ethanol, and 400 mL of distilled water.The preparation was carried out in glass bottles of 1 L capacity.The final mixture was sterilized in an autoclave at 121 °C and 15 PSI for 20 min.Using a laminar flow hood, the sterile medium was poured into 50 mL Falcon ™ polypropylene tubes until completing 15 mL, allowing it to dry for 12 h.Finally, the tubes were kept at 26 °C and 60% relative humidity until the day of the experiment.

Experimental rearing conditions
In Scolytinae with haplodiploid reproduction, adults can mate with siblings within the natal nest before dispersal.Based on this biological characteristic, we assume that all the females used for the experiment were fertilized 130,175,176 .Prior to inoculation, X. bispinatus females were sterilized by immersion for 5 s in 70% alcohol, and later, in water for the same time to remove the alcohol.Surviving females were individually inoculated into 50 mL Falcon™ tubes containing culture medium.To favor gas exchange and prevent individuals from escaping, the tube caps were perforated and covered with a metal mesh.The inoculated tubes were placed in rearing chambers at constant temperatures of 17, 20, 26, 29, and 35 °C, arranged vertically and kept in complete darkness (Supplementary Fig. S1 online).Each temperature treatment included 90 tubes, resulting in a total of 450 inoculated tubes for the entire experiment.Temperature monitoring in each growth chamber was performed using the Elitech URC-4 device.

Experimental count of individuals
Every four days (throughout 36 days that the experiment lasted), 10 tubes from each rearing chamber temperature were checked.The artificial medium was carefully dissected (destructive monitoring) and the number of eggs,

Egg hatching percentage
Accounting for hatched eggs is extremely important when estimating the growth rate of the species.However, due to the cryptic habit of X. bispinatus, which requires destructive monitoring of the culture medium, it was not possible to make such an estimate.To correct this, we used information from Dendroctonus ponderosae, a species related to the genus Xyleborus 79,102,177 .We consider that this procedure represents a more viable alternative than assuming a 100% hatching rate.The use of proxies or assumptions to obtain limited demographic information is a valid method, especially when it is necessary to make management decisions before they can be collected 178,179 .Furthermore, it is important to highlight that Xyleborus bispinatus present similar patterns regarding optimal growth and the upper and lower population thresholds (Fadda et al., in review) as D. ponderosae when raised at similar laboratory temperatures 180,181 , thus enhancing the robustness of our methodological approach.
The information regarding the percentages of hatched eggs was obtained from mountain pine beetles (Dendroctonus ponderosae) reared axenically on a standardized diet at constant temperatures of 10, 18, 20, 24, 27, 32, and 35 °C180 .Using this information, we fitted a second-order polynomial curve using the nls.lm function of the minpack.lmR package 182,183 .From this curve, we obtained an estimation of the egg hatching percentage equivalent to the temperatures tested in our work.Subsequently, these values were used in conjunction with the number of egg laying made by X. bispinatus to estimate the net reproductive rate necessary to calculate the growth rate of the species at its different temperatures (Eq.1; Supplementary Fig. S3 online).
where p = Proportion of hatched eggs, T = Temperature, a, b and c = parameters of the regression.

Survival of developmental stages
Survival ( S k ; Eq. 2) was calculated for each of the stages of the life cycle of X. bispinatus (see section "Experimental count of individuals").This is computed by the quotient of the individuals belonging to the immediate following stage of development to the one of interest (k + 1) and the value obtained for the stage that we wish to estimate (k).The equation returns values between zero and one, and is estimated as follows: where S k = Survival at development stage k, N k+1 = Number of individuals at stage k + 1 , N k = Number of individuals at development stage of interest.
In the case of eggs, the hatching percentage at each temperature was first estimated (Eq. 1) and then the survival rate.In adults, survival was only estimated for females since they are mainly responsible for the population growth of the species 67 .The number of individuals was estimated from the quotient between the total number of females (live and dead) and the values of live females recorded in the experiment.

Net reproductive rate
The net reproductive rate ( R 0 ; Eq. 3) was estimated according to what was stipulated by Southwood 184 as the product of survival in each stage of development and the proportion of hatched eggs estimated in section "Egg hatching percentage" at each temperature: where R 0 = Net reproductive rate, p = Proportion of hatched eggs, H = Number of eggs counted, S k = Survival rate at stage k.

Population intrinsic growth rate
The intrinsic growth rate of the population ( r ; Eq. 4) was estimated as the quotient between the natural logarithm of the net reproductive rate and the generation time ( G ) observed at the temperature of interest.G is defined as the time between the first observation of the egg stage and the first count of an adult hatchling of the progeny in question at a given temperature: where (r) = Population intrinsic growth rate, R 0 = Net reproductive rate at a given temperature, G = Generation time, days from the first observation of eggs to the first record of females from the progeny.
(1) We fitted a non-linear model to the values of the intrinsic growth rate obtained for each temperature (section "Population intrinsic growth rate"), using the nls.lm function of the minpack.lmR package.This was based on a convex mathematical function (Eq.5) whose parameters are biologically meaningful: an optimum suitability for the species (centroid), expressed in terms of optimum temperature t opt where the maximum intrinsic growth rate of X. bispinatus occurs (r max ) , and the kurtosis of the fitted curve ( v o ; a positive real number that determines the niche breadth): r = Intrinsic growth rate, r max = Maximum intrinsic growth rate, T = Temperature, T opt = Optimum growth temperature, V 0 = Kurtosis of fitted curve.

Ecological niche modeling
To mitigate the fact that we modeled the intrinsic growth rate as a function of one variable (temperature) and this parameter might depend on other factors such as humidity 32 , we applied a two-step modeling approach before evaluating the intrinsic growth rate in the geographic space.First, we employed correlative niche models fitted using bioclimatic variables to consider the effect of seasonal variations of precipitation and temperature on the potential distribution of this species; the above allowed us to delimit the sites where the evaluation of r made biological sense 185 .Then, we evaluated the function that relates the mean temperature with the intrinsic growth rate (Eq.5) on those sites with suitable conditions of precipitation and temperature.

Occurrences
We compiled georeferenced and dated records of X. bispinatus in America from: BarkBeetles (https:// www.barkb eetles.info/), the Global Biodiversity Information Facility (GBIF; https:// www.gbif.org/) 186 , bibliographic review of scientific journals, and collected specimens at Los Tuxtlas, Veracruz.With Google Earth Pro 7.3, we verified these records coincided with field collections and discarded records from cities, museums, and other scientific collections.Additionally, we eliminated temporal (same year) and spatial (at 30 arc seconds) duplicate records.This produced a total of 106 presence points for the species of which 70% were used for calibrating the model and 30% for its validation from a random partition.

Environmental variables and time specific niche modeling
Monthly climatic data on minimum temperature, maximum temperature and cumulative precipitation were obtained from the CHELSAcruts database 187,188 at a resolution of 30 arc seconds for the period comprising 1986-2016.With this information and following the methodology in O'Donnell and Ignizio 189 , bioclimatic layers were constructed for each year during the period in question.
To estimate the ecological niche of X. bispinatus, we used time-specific niche models (also known as time calibrated species distribution models 190 ).This approach allows us to calibrate niche models with spatial information of high temporal resolution, reducing niche and distribution estimation biases 191 .The modeling algorithm we employed was minimum volume ellipsoids, as they are a simple and biologically informative way to represent an n-dimensional hypervolume according to the original proposal of Hutchinson 11 .This proposal considers a niche to have a convex shape where the maximum fitness is at the centroid of the ellipsoid and decreases towards the periphery 16 .Moreover, it has been documented that the distance to the center of an ellipsoid is related to fitness attributes, such as population abundance and genetic diversity 192,193 .
The modeling process was carried out using the tenm (temporal ecological niche models) 194 package in R 183 , which allows a time-specific model selection process to be executed.In this process, the data were first thinned using a distance of 0.0083 degrees to reduce problems associated with spatial autocorrelation.Likewise, the environmental information corresponding to the year of the presence records was extracted and the Pearson correlation coefficient was estimated; those whose correlation was less than 0.8 were chosen as modeling variables.The selection of better models was done considering a significance of 95% for the partial ROC statistical test 195 (1000 iterations and 50% of the data was used for the bootstrap), in addition, the partial AUC ratio and AUC value were calculated from 50,000 randomly selected environment points, a performance criterion of 10% omission and with the lowest environment prevalence percentage to decrease the potential error of commission and thus avoid overestimation; this selection process has been used in other modeling works 192,196,197 .Finally, these were converted to a binary prediction (presence-absence) based on a 10th percentile threshold.

Adjusted growth map from the intrinsic growth rate function
Prior to running the models, we designed as the calibration area 185 the geographical region encompassing the Mexican territory together with the sympatry zone shared by the species X. glabratus and X. bispinatus in the United States, given that the first might be infected by the last.Then, in order to perform the metapopulation process-explicit dynamic models, we estimated the growth rate in the study area by evaluating Eq. 5 on each pixel with the annual mean temperature layer (2.5' resolution; ~ 5 km per pixel) from WorldClim 120 ; this produced a r i map.To consider the effect of other environmental variables such as precipitation on growth rate and to reduce the overprediction related to some sites having temperature values outside the fundamental niche of X. bispinatus at certain times of the year, the r i map was multiplied by the binary suitability map of the final model.Finally, with the crop and mask functions of the raster R package 198 , we proceeded to delimit the polygon where the simulation was carried out.www.nature.com/scientificreports/

Modeling the invasion dynamics of Xyleborus bispinatus
The spatiotemporal dynamics of X. bispinatus was estimated as a function of climatic conditions, growth rate and density dependence factors.For this, the metapopulation model of Eq. 6 was used.In this model it was assumed that the species grows in a region divided by a regular grid with con i = 1, 2 . . .n number of cells.In the absence of migratory processes, the growth of a population x i is determined by its biology and the environmental characteristics of cell i .During the invasion process, the populations x i are connected by dispersal, therefore, the population abundance in each cell i also depends on the rates of entry and exit of individuals in cells i .The metapopulation model is shown below 16 where ẋi is the rate of population change in cell i at time t .r i is the intrinsic growth rate of cell i which was obtained from the growth rate map for the zone of interest.a is a dense-dependence factor which measures the intensity of intraspecific competition.δ hi is the immigration rate from cells h to cell i ; δ ih is the migration rate from cell i to the other cells.The immigration and migration rates are determined by a distance-dependent first- neighbor dispersal kernel (Eqs.7A and 7B).This procedure composes a quadratic system of coupled nonlinear differential equations that requires solving as many equations as the number of cells, which implies a relevant processing time as well as a large computational capacity (deciphering approximately 63,000 equations at the resolution and spatial extent worked).A resolution of 2.5 min per cell was considered for the simulation in the calibration area, since at finer resolutions the numerical solution of the model was not possible.The value of the density dependence parameter a was set at a = 0.0000001, to produce population densities (in units of individuals per km 2 ) of approximately 130,000 individuals in the most suitable cells (considering that the carrying capacity in the model without dispersion is r max a ).

Dispersion kernel
To model dispersal between cells, we used an exponential kernel 50,199 where D max is the maximum distance an individual can travel per day between patches, and was established at ~ 5 km 67,162,163,200,201 ; δ ih and δ hi are migration rates of individuals between cells; k is the maximum migration capacity, which was set as one tenth of the intrinsic growth rate r i and b is a positive constant that modulates travel capacity.The distance between cells was measured with Euclidean distance using the cell coordinates h and i .The barriers that prevent migration from patch h to patch i are simulated by taking δ hi = δ ih = 0 .In our case we relied on an irreducible dispersal matrix; that is, www.nature.com/scientificreports/there is the possibility of reaching all grid elements, from any other cell including those of low suitability and without taking into consideration biotic barriers such as the Allee effect 19 .

Starting points for modeling invasion routes
We defined various starting points in Mexico to conduct the simulations, including ports and customs.This decision was primarily influence by the fact that the great majority of ambrosia beetles are favored by anthropogenic dispersal, which facilitates their escape from natural barriers that might otherwise confine them to their native ranges.Thus, we considered the geographic coordinates for the ports of Salina Cruz, Veracruz, Manzanillo, Altamira and Ensenada, and the land customs of Nuevo Laredo and the Nogales.In addition, we also considered a site in Texas where X. bispinatus and X. glabratus are geographically close and the disease is known to be present, representing a higher risk for Tamaulipas as a potential gateway through which the vector and the pathogen could enter Mexico by land (Supplementary Table S1 online).The summary diagram of the workflow carried out in this study is shown in Fig. 6.

Ethics approval and consent to participate
This article does not describe any studies involving human participants performed by the authors.All applicable international, national and/or institutional guidelines for the care and use of animals were followed.

Figure 1 .
Figure 1.Mathematical function of growth rate of X. bispinatus at different temperatures considering a 95% confidence interval.

Figure 2 .Figure 3 .Figure 4 .
Figure 2. (A) Map with the suitable sites estimated from the ecological niche modeling, the records of X. bispinatus presence, and the starting coordinates used for the spatiotemporal simulations.(B) Projection in geography of X. bispinatus growth estimated from the intrinsic growth rate function (r i map) used for the construction of the process-explicit models.This figure was generated in ArcGIS ver.10.4 (https:// www.arcgis.com/ featu res/ index.html).

( 6 )Figure 6 .
Figure 6.Workflow diagram used for the construction of the process-explicit model simulations.The figure was assembled with the growth curve depicted in Fig. 1, along with maps created using the raster package version 3.6 in R and ArcGIS version 10.4.(https:// www.arcgis.com/ featu res/ index.html).
, pupae and adults (males and females; Supplementary Fig.S2online) was counted.In addition, for this last stage of development, living and dead individuals were recorded.Observations were made with a Leica EZ4 stereomicroscope.The experiment took place over 36 days, totaling nine colony monitoring openings.Individual counts at each colony opening date were averaged and accumulated.With this final accumulated average of the experiment for each stage of development and temperature, the biological parameters of the species were estimated.